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Abstract 

A successive continuation method for locating connecting orbits 
in parametrized systems of autonomous ODEs is considered. A local 
convergence analysis is presented and several illustrative numerical 
examples are given. 
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1 Introduction 

The existence of a trajectory connecting equilibrium points of an ODE, a 
homoclinic orbit or heteroclinic orbit, also called connecting orbit is of signif- 
icance in a variety of applications. Connecting orbits often arise as traveling 
wave solutions of parabolic and hyperbolic PDEs., e.g., in combustions mod- 
els p|. They have been shown to underlie intermittency phenomena in fluid 
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mechanics [H , "bursting" in models of biological cells |BU| , chaotic vibration 



of structures [29], and chaotic behavior of electronic circuits, P, [181 19], light 



pulses in fiber optics ||28|| , and chemical reactions [p4}[ . The corresponding 
numerical problem is that of finding solutions (u(t), A) of the system of au- 
tonomous ODEs 

u'(t) - f{u(t), A) = 0, <), /(•, ■) G W 1 , A G R»\ (1.1a) 
lim «(/!;) = w , lim w(t) = u\. (1.1b) 

t— »— oo t— >+oo 

Most algorithms for computing connecting orbits reduce (1.1) to a boundary 
value problem on a finite interval using linear or higher order approximations 
of stable and unstable manifolds near u and ui, respectively. We studied and 



applied such a method in [13[ and generalized it in ||21|| . An alternate conver- 



gence analysis of the method is given by Schecter [R3] . The basic method of 



choosing appropriate boundary conditions dates back to Lentini and Keller 
J27]. A closely related procedure, which uses generalized eigenspaces for 
fu(u , A) and A) to construct the relevant projections, was developed 

by Beyn P, || and used by Champneys and Kuznetsov in their recent work 
on locating higher codimension homoclinic orbit bifurcations 171] . A shooting 
method for computing connecting orbits has been used by Rodriguez-Luis 
et al [3T|. For parameter superconvergence results see the recent work of 



Sandstede 32 



We considered the case of center manifold in |Z0], . Champneys, Kuznetsov, 
and Sandstede have also generalized the algorithm in [7j to the case of non- 
hyperbolic equilibria ||. Both, hyperbolic and nonhyperbolic equilibria are 
considered by Canale ||, using multiple shooting for discretizing the approx- 
imate truncated problem, together with a new a new technique for evaluating 
partial derivatives in the boundary conditions. 

In this paper we present a local convergence analysis of the successive 
continuation method for locating connecting orbits. This method aims at 
computing a connecting orbit when there is a region of uncertainty for the 
parameter values, which is the case in most practical problems. In contrast, 
most previous work is concerned with the continuation of connecting orbits, 
given a sufficiently accurate starting connecting orbit. 

Assume that uq and u\ are hyperbolic fixed points of the vector field /, 
and that the local unstable manifold W£ c (uo) at uq has dimension no, and 
that the local stable manifold Wf oc (ui) at u± has dimension ri\. Let So and 
Si denote the tangent planes to W^Juq) and Wf oc (ui), respectively. In the 
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finite interval approximation, the scaled independent variable takes values 
in [0,1]. A general outline of the full algorithm for locating and continuing 
connecting orbits is then as follows: 

(I) Time integration to get an initial orbit. Typically, this gives a very 
inaccurate "approximation" of a connecting orbit with either (u(0)—uq) G Sq 
or (u(l)—U\) G Si. This step can be done by continuation with the "period", 
i.e., the total integration time, as parameter. 

(II) Locate a connecting orbit. Use a sequence of homotopies ("successive 
continuation") that lead to an orbit with both (u(0) — uq) G So and (u(l) — 
Mi) G Si. 

(III) Increase the accuracy. Use continuation to further increase the pe- 
riod of the approximate connecting orbit. 

(IV) Continue the connecting orbit. If desired, compute an entire family 
of approximate connecting orbits. 

In this paper we are primarily concerned with Step (II). In Section 2 we 
formulate the approximate problem and a successive continuation algorithm. 
This algorithm is implemented in a general code based on the continuation 
package AUTO |17| and is used to compute the examples in Section 4 below. 
In Section 3 we prove a local result, namely that the successive continuation 
algorithm has an open neighborhood of convergence. More precisely, we show 
that for any point in this neighborhood there exists a continuous piecewise 
smooth homotopy leading to the solution. Section 4 contains examples. In 
Section 5 we conclude with a discussion of the applicability and limitations 
of the algorithm. 

In our original algorithm for continuing connecting orbits [[il| |2(| we 
included the eigenvectors as unknowns in the full system to be solved at each 
continuation step. This can lead to an ill-conditioned system in the case of a 
nearly defective eigenvalue, eigenvectors becoming nearly linearly dependent. 
In the present paper we use orthonormal bases to construct the appropriate 
projections associated with the real Schur factorization [2J] of f u {uo, A) and 
f u (ui, A) or their transposes. This procedure is more accurate and stable. 

We first realized the power of the successive continuation approach for lo- 
cating connecting orbits when trying to find a homoclinic orbit in a singular 
perturbation problem |T5|| . For application to a singularly perturbed sine- 
Gordon model of the Josephson junctions and to the Ho dgkin- Huxley equa- 



tions see ||14|| . The basic idea was also used by Champneys and Kuznetsov 



to locate homoclinic orbit in Chua's electronic circuit and in the FitzHugh- 
Nagumo equations J7|. In the context of optimization, the approach is also 
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used in the tutorial papers |TH . 



2 Description of the algorithm 

In addition to the requirement that the fixed points Uq and u\ be hyperbolic, 
we assume that the eigenvalues of f u (u , A) and f u (ui, A), respectively, satisfy 

Re/i 0n < ... < Re/i 0no+1 < < fi 0A < Re/i 02 < ... < Re/i 0no ,(2.1a) 
Re/x xl < ... < Re/x lni < < Re/i lni+1 < ... < Re/i ln (2.1b) 

The method can be extended to the case fi 01 = 0, as in || It also 
extends to the cases of complex and multiple fi 1 by a modification of Step 0, 
Eq. (2.13) of the algorithm (see Section 4.3 for an example). The algorithm 
requires evaluation of various projections onto tangent subspaces at uq and 
Mi. We construct these projections using the real Schur factorizations |23 



fu(u ,X) = QoT Qq, (2.2) 

and 

fu(u 1 ,\)=Q 1 T 1 Ql (2.3) 

of f u (u , A) and f u (ui, A), respectively. The factorization (2.2) has been cho- 
sen so that, in the real case, the eigenvalues fi 01 , /x no appear in the upper 
left corner of T . In the complex case the corresponding 2x2 blocks for 
each complex congugate pair appear. A similar choice applies to the eigen- 
values 1; ...,/i lrtl and T\ in the factorization (2.3). Qo = [<7o,i--9o,ra] an d Qi 
— [9i,i"-9i,n] are orthogonal, and T and T\ are upper quasi-triangular (1-by-l 
and 2-by-2 blocks on their diagonals). Then the first n columns g ,i; Qo,n 

of 

Q form an orthonormal basis of the right invariant subspace S of f u (u , A), 
corresponding to /i 01 , fi 0)Tl0 , and the last n — n columns go,n +i; Qo,n of 
Qo form an orthonormal basis of the orthogonal complement S$. Similarly, 
the first rii columns qi i, qi tU of Qi form an orthonormal basis of the right 
invariant subspace Si of f u (ui,X), corresponding to 1; /x 1 ni , and the 
last n — rii columns ?i, m +i, •••,?i, n of Qi form an orthonormal basis of the 
orthogonal complement S^. 
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Approximate problem. 

The approximate finite interval problem is now as follows: given eo = 
e* G R+, "small", T G R+, "large", find a solution 

(u , A , -u , d , d ± , e 1 , Q , Q 1 , T , T x ), 

where u* G C^fO, l]),R n ), A* G R" A , ujj, m*, d^, d\, G M", e* small, Q* ,Ql,T*,T* 
W l x ]R n , of the time-scaled differential equation 

u'(t)-Tf(u(t),\)=0,0<t<l, (2.4) 

subject to stationary state conditions 

/(«o,A) = 0, (2.5a) 

f( Ul ,X) = 0, (2.5b) 

left boundary conditions 

u(0) = u + e d , (2.6a) 

I d |= 1, (2.6b) 

d ■ qo,n +j = 0, j = 1, n - n , (2.6c) 

right boundary conditions 

= Mi + ei^i, (2.7a) 

I d 1 |= 1, (2.7b) 

Tj = dj -gi.m+j = 0, j = l,...,ra-ni, (2.8) 

with real Schur factorization of f u (u , A) 

/ u (u , A)g 0) < = 9oj'*S,i' * = !> ( 2 - 9a ) 

conditions for Q to be orthogonal 

9ji9oj = ^j, i = 1, j, j = 1, ...,n, (2.9b) 
real Schur factorization of f u (ui, A) 

n 

f u (u l7 X)q lti = 9ij*j,i> i = 1, -, n, (2.10a) 



and conditions for Qi to be orthogonal 



9i,i9iJ = s ij> 1 = ij-jJj i = (2.10b) 

Here the matrix T$ = [tj J, 5 = 0, 1, is upper quasi-triangular with all ele- 
ments below the subdiagonal equal to zero. During each continuation step 
its n(n + l)/2 elements above the subdiagonal vary, while the remaining n — 1 
elements on the subdiagonal are either fixed at zero or vary. In the latter case 
nearby diagonal elements are equated. Specifically, we have the following: 
each 1-by-l block t\ i on the diagonal is a real eigenvalue, and each 2-by-2 



block on the diagonal 



t 5 t 5 

b i+l,i 



corresponds to a complex conjugate 



pair, where tf +l i+1 = tf { is the real part of the eigenvalue and tf i+1 = —tf +lji 
is the imaginary part. In the process of continuation we have the following 
special cases: (i) if all eigenvalues are real and distinct then we let all entries 
on and above the main diagonal vary, while the entries on the subdiagonal 
are fixed at value zero; (ii) if we have a complex conjugate pair then we add 
the equation t 5 i+l i+1 = t\ { and let tf +l i vary, (iii) if we have a double real 
eigenvalue then we can choose either of the preceding options. To summarize, 
we have 

if t° i+hi+1 ± t% set tl ht = 0, else set f° +1>i+1 = t% and vary t° i+1>i , 

i = l,...,n-l, (2.9c) 

and similarly 

if + t],i set t] +hi = 0, else set t\ +w = t]- and vary t} +1>i , 

% = l,...,n-l. (2.10c) 

We remark that a drawback in the earlier version of our algorithm was its 
inability to compute through double eigenvalues: we had to stop at double 
real eigenvalues and restart with a newly -defined continuation algorithm. 
The present method overcomes this problem: one just fixes appropriate pa- 
rameters and frees other appropriate parameters (cf. steps 3 and 4 in Section 
4.2). 

The above constitutes n differential equations with n c = 3n 2 + 7n + 2 — 
no — n\ constraints and n v = 3n 2 + 5n + ri\ + 1 scalar variables. Generically, 
it is necessary that n c — n v = n, in order for the system (2.4) - (2.10) to have 
a unique solution. This gives the relation n\ — n — (n + ni) + 1. 
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By a slight modification of the analysis in JH]], using a version of the 
implicit function theorem used, e.g., in Beyn [||, it is easy to show, un- 
der appropriate transversality conditions, for T large enough and eg small 
enough, that the system (2.4) - (2.10) has a unique solution which is a good 
approximation to the solution of (1.1). Note that since in the present case we 
look for only one connecting orbit rather than a branch of connecting orbits, 



the dimension of the parameter vector A is one less than in [21|. We also note 



that the proof in |HJ] is for the approximate problem in the same Banach 



space as the exact problem and that the result is independent of the partic- 
ular implementation, and hence applies to the algorithmic implementation 



(2.1) - (2.6b) in [21] as well as to (2.4) - (2.10) in the present paper. 
Algorithm for locating a connecting orbit. 

The solution to the above system is found via a sequence of homotopies 
that locate successive zero intercepts of the Tj, in 

Tj - d-t ■ q ljni+ j(u!, A) = 0, j = 1, n - n x , (2.11) 

(cf. equation (2.8)). In each homotopy step we compute a branch, i.e., a one- 
dimensional manifold, of solutions. For this we must have n c — n v = n — 1, 
and hence n\ = n — (no + n\) + 2; ri\ > 0. 

Let So.fcj k — 1, no, be the right invariant subspace of f(uo, Xq) corre- 
sponding to the eigenvalues fj JQ1 , fi ok . Then the first k columns go,i> Qo,k 
of Qo form an orthonormal basis of Sq&. Initially we replace (2.6) by the 
equivalent equations 

no 

u(0) =u + e ^2 c jQoj, (2.12a) 

no 

E c ? = 1 - ( 2 - 12b ) 

Step 0. Initialize the problem parameter vector A, and set the algorithm 
parameters e and T to small, positive values, so that u(t) is approximately 
constant on [0, T]. Set do = goi? 

u (t) = u + e d , < t < 1, (2.13) 

ei = \u(l) — U\ \ , d\ — (u(l)— «i)/ei, C\ = 1, and c 2 = ... = c no = 0. Here the 
initial direction of the orbit is typically chosen along the eigenvector of the 
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weakest unstable eigenvalue, d = g i? which implies (2.6), i.e. u(0)— u G 5*0. 
Step 1. Compute a solution branch to the system (2.4), (2.7), (2.11), 
(2.12a), in the direction of increasing T, until u(l) reaches an ei-neighborhood 
of ui, for some e\ > 0. Scalar variables are T, e\ G R, di 6 1", r £ R" - ™ 1 . 
There are n differential equations with n c = 3n — n\ + 1 constraints and n B = 
2n — ni + 2 scalar variables, and hence n c — n v = n — 1. This initial solution 
is normally a very inaccurate approximation of the solution of (2.4) - (2.10) 
since, in general, the Tj in (2.11) will be nonzero and hence u(l) — u\ ^ S\. 
In practice one typically continues until e± stops decreasing. Its value is then 
not necessarily small, but the successive continuation procedure is intended 
to work even then (see also Section 5). 

Step 2 (for n > 1 ), k = 2. Keep T fixed and free C\ and c 2 , i.e., we let 
u(0), which had been confined to So,i in Step 1, now run through So,2- More 
precisely, we compute a branch of solutions to the system (2.4), (2.7), (2.11), 
(2.12) to locate a zero of, say, t\. Free scalar variables are e± : ci, c 2 G R, di G 
R n , t G R n ~ ni . There are n differential equations with n c = 3n — rii + 2 
constraints and n v = 2n — rii + 3 scalar variables, and hence n c — n v — n — 1. 

(for n > 2 ), k = 3, n . Keep T and Tj = 0, j = 1, ...k — 2, 
fixed. Also fix one more component of r, say Tfc_2, which reached value zero 
in the previous step. Thus we now let -u(O) — uo, which had been confined to 
«So,fc-i, run through So^- More precisely, we compute a branch of solutions 
to the system (2.4), (2.7), (2.11), (2.12) to locate zero of, say, r fe _i. Free 
scalar variables are ei, c±, c&, t^-i, T„_ m G R, d\ G R n . There are n 
differential equations with n c = 3n — n\ + 2 constraints and n B = 2n — ni + 3 
scalar variables, and hence n c — n v — n — 1. 

In the following steps, A also varies. For the purpose of implementation, 
(2.12) is better than (2.6), but the latter is more stable when eigenvalues 
coalesce. There are examples where continuous dependence of goj on A is 
lost, and where Cj oscillates. Consequently, starting with Step 3, we use 
(2.6) instead of (2.12), since (2.6) is more stable in this case. 
Step 3, k = n +l, n +riA= n — ni+1. The free parameters ci, c*, G R 
are replaced by d G R n . Fix another component of r, say, r fe _ 2 , which reached 
zero in the previous step, and free one more component of A, say, Xk- no - More 
precisely, we compute a branch of solutions to the system (2.4) - (2.7), (2.9) 
- (2.11). There are n differential equations with n c = 3n 2 + 7n + 2 — n — rii 
constraints and n v = 3n 2 + 6n + 3 — no — ri\ scalar variables, and hence 
n c — n v = n — 1 . 

In some of the examples below we use slight variations on the basic algo- 
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rithm in order to minimize the number of steps, which attests to its flexibility. 

The above algorithm corresponds to Step (II) in the outline of the full 
algorithm for locating and continuing connecting orbits in Section 1. Step 
(IV) of that outline can now be described as follows. 

Algorithm for continuing a connecting orbit. 

Compute a branch of solutions to the system (2.4) - (2.10), where T, and 
all components of A vary. Here w(0) G 5*0, u(l) G Si, T, and all components 
of A vary. A phase condition 

[\u'(t)-q'(t)) ■ u(t) dt = (2.14) 
Jo 

may be added if T is kept fixed and e and e\ are allowed to vary. Here q(t) 
is a previously computed orbit on the branch. 

Remark 1 The factorizations (2.2), (2.3) can be computed at each contin- 
uation step using, for example, LAPACK, or by continuation. Continuation 
appears to be more robust. This is to be expected when the eigenvectors, and 
hence the Schur matrix, vary rapidly with a slow change of a continuation 
parameter, since in this case the pseudo-arclength step will automatically de- 
crease to capture the transition. We also compared computing a connecting 
orbit via continuation and via factorizations in cases where the Schur ma- 
trices depend "smoothly" on the problem parameters. Computations using 
continuation were 10-20 times faster (due to larger continuation steps). 



Remark 2 For the convenience of the reader, we summarize here some of 
the issues related to the algorithms for locating and continuing a connecting 
orbit, in this paper and in our earlier work. (1) The Approximate prob- 
lem in here is to locate a connecting orbit, while the approximate problem 
in [ 2l J is to continue a branch of connecting orbits. Hence the difference be- 
tween two formulations. (2) Once a connecting orbit has been located, there 
are two basic methods to continue a branch of connecting orbits (see Algo- 
rithm for continuing a connecting orbit on p. 7): (a) T, and A vary, 
while e and €\ are kept fixed; (b) T is kept fixed, eo and e\ are allowed to 
vary, and a phase condition is added. Both types of continuation have their 
advantages and disadvantages. The advantage of (a) is that the accuracy of 
the computed orbits on the branch is the same (since eo and €\ are fixed), 
while the advantage of (b) is that it often works (and works well) for more 
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difficult cases (i.e. when the orbit changes rapidly during the continuation), 
when (a) fails or becomes extremely slow. See also the Remark in [T3jJ which 
explains (b) in more detail. 



Remark 3 An alternate formulation. An alternate formulation using 
transposed matrices, as in ^ 0, 0, W\h results in a similar algorithm. The 
real Schur factorizations (2.2) and (2.3) are replaced by 

f*{u Q ,\) = Q Q T Ql (2.15) 

and 

ff(u 1 ,\) = Q 1 T 1 QT, (2.16) 

respectively. Correspondingly, the system (2.4) - (2.10) is replaced by the 
system (2.4) - (2.6a,b), (2.7), and 

do ■ qo,j = 0, j = 1, ...,n-n , (2.17) 

Tj = di ■ qi tj = 0, j = 1, ...,n- m, (2.18) 

n 

fu( u o,X)Qo,i = ^2%,j t0 j,v i = l,...,n, (2.19a) 
i=i 

*/ ^ t°i then set t°i+i,i = °> else set = and var V *°+i,i> 
% = 1, ...,n-l, 

9ji9o,j = « = 1, j = l,...,n, (2.19c) 

n—ni 

/J(«i,A)gi,i = ^gij^, i = l,...,n, f2.20aj 
i=i 

i = l,...,n-l, f&gfl^ 

= ^J' z = l,-..,J, J = l,.»,n. (2. ,20c j 
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3 Convergence of the algorithm 



The following notation will be used throughout this section: 

ri\ = n + l — n — ni, (3.1a) 
n T = n — n\ — n\ + no — 1. (3.1b) 

In the system (2.4)-(2.10), equations (2.9)-(2.10) can be used to find nu- 
merically Qo and Q 1 as functions of u , A and tii, A, respectively. This allows 
to treat (2.4)-(2.8) as a self-contained problem. Replace the problem (2.4)- 
(2.8) by an equivalent one 

u'(t)-Tf(u(t),X)=0, 0<t<l, (3.2) 



f(uo, A) = 0, (3.3a) 
f(ui, A) = 0, (3.3b) 

|«(0) - u | = e , (3.4) 

Oi = (u(0) - u ) • qo, no +i(uo, A) = 0, i = 1, ...,n - n , (3.5a) 
o"n-n +i = (u(0)-u )-q ,i+i(u ,\)-iy i = 0, i = 1, ...,n - l(3.5b) 

r i = — u i) ' 9i,m+i( u i) A) = 0, i = l,...,n T . (3.6) 

The equivalence is apparent, since the system (3.2)-(3.4), (3.5a), (3.6) be- 
comes formally equivalent to the system (2.4)-(2.8) if appended by the equa- 
tions 

d = (u(0) - u ) /e , (3.7a) 
ei = \u(l) — tii| , (3.7b) 
di = (ti(l) - tii)/ei, (3.7c) 

whereas (3.5b) merely defines the parameters i/j. Notice that Ui = Cj+i, i = 
1, ...,rio — 1, where q's appear in (2.12). 

Definition A set A in a Banach space X will be called an arc if there 
exists a Banach space Y and an open sets U C X, V C Y and a C 1 - 
diffeomorphism F : U — > V such that A C U and F(A) is a closed bounded 
rectilinear segment (a degenerate segment consisting of a single point is al- 
lowed). 
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Remark 4 Any arc A has finite length. Indeed, if 7 : [0, 1] — > Y is a 
linear parametrization of the rectilinear segment F(A), then F~ x o 7 is a 
smooth parametrization of A with continuous, hence bounded, velocity v(s) = 
[F" 1 ] ( 7 (s)) o y( s ). Thus (length of A) < sup \v(s) \ . 

Introduce two Banach spaces, X = C^O, 1] x R n x IR n x R Ut = {x = 
(U,A) : U = (u,uo,U!) G C^l] xR n xW l , ue C^O, l],u ,«i G M n , A = 
(1/, A) G M n 7 1/ G IT 0-1 , A G M nA } and Y = C[0, 1]xR"xK"xRx M n_1 . 
Define F : X -> Y by 

F([/,A) = (u'-T/(u,A), f(u ,X), f(u±, A), |u(0)-u | 2 -ejj, <r(l/,A)), 

(3.8) 

where a : X — > M n 1 , a = (a"i, a" n -i) is defined by 
o-i(U,A) = (u(0) - u ) ■ qo, no+ i(u ,X), i = l,...,ra-ra 

(3.9a) 

(Ji(C/,A) = (u(0) - W ) • g ,i+l-(n-no)(«O> A ) - ^-(n-no), 1 = W - Wo + 1, W - 1. 

(3.9b) 

Also define r : X -> R% r = (n, r n J by 

T j{U,A) — (u(l) — tii) • A), j = l,...,n T . (3.10) 

Theorem Dearie F^ :I^7x IR nT , F^ = (f^\f^\ ...,F^ , j = 
l,...,n T , 6y 

F^([/,A) = (F^A)^!^^),..,^^)^^!,..,^), j = l,..,n T -l, 

(3.11a) 

FM(U,A) = (F(U,A),r(U,A)). 

(3.11b) 

Assume that there exist (U*,A*) G X s^c/j t/iat 

F (nT) ([T,A*) = 0. (3.12) 

Assume also that for j = l,...,n T , F^ is a ^-mapping in an open neigh- 
borhood of (U*,A*), and its Frechet derivative at (U*,A*) is an invertible 
bounded linear operator. 



12 



Then there exists an open neighborhood V of (U*,A*) such that for any 
(C/°,A°) G V H {([/, A) : F(U,A) = 0} inere exists a continuous piecewise 
smooth curve x(s) = (U(s),A(s)), s G [0,1], x(0) = (C/°,A°), x(l) = 
(£/*, A*) u^/j t/ie following properties: 

(i) i([0,l])cV and consists of n T arcs, specifically, there exist numbers 
< s < si < .. < s Ht = 1 such that each „4 fe = x([sfc_i, s fc ]), k — 1, n r is 
an arc; 

/or any s G [0, 1], F(x(s)) = ; 
(mj for any k G {2, ...,n T }, if x G »4fc, #ien = ... = Tfe_i(x) = 0; 

/or any s G [0, 1], if Tk(x(s)) ^ /or some k G {1, n T — 1}, taen 
A fc+ i(s) =A° +1 ,...,A nr (s) =A° t . 

The hypotheses of the Theorem which require that Frechet derivatives of 
F^ at the solution (U*,A*) are invertible bounded linear operators, express 
the following transversality requirements. The surface {F = 0} and the 
hypersurfaces {t\ = 0},...,{r„ T = 0} intersect transversely at the solution 
(U*,A*); the same holds for the surface {F = 0} and the hypersurfaces 
{n = O},..,-^ = 0},{A j+1 = A* +1 },-.,{A nT = A;j, for each j = 1, ...,n T - 
1. In particular, these requirements imply that, at the end of each step j 
(j — 1, ...,n T ), Tj crosses zero transversely. 
Define 

i 

E = {iGl: F(x) = 0}, Si = S n f^\{x G X : T k (x) = 0}, i = l,...,n T . 

k=l 

(3.13) 

The following Lemma 1 may be interpreted as Lemma 2j below for j = n T . 
To emphasize this, the condition (iv) is included in Lemma 1 even though 
it is a trivial consequence of (iii). Incorporation of the case j = n T into the 
statement of Lemma 2j would lead to a cumbersome formulation. 

Lemma 1 Let F, r, U* , A*, F^ satisfy the hypothesis of the Theorem. 
Then there exists an open neighborhood W of (U*,A*) such that for any 
(U nT ~ l , A nT_1 ) G W^nS„ T _i there exists a smooth curve x(s) = (U(s),A(s)), 
s G [0,1], x(0) = {U riT ~ 1 , A nT ~ r ), x(l) = (U*,A*) with the following proper- 
ties: 

(i) x([0, 1]) C W and is an arc; 

(ii) for any s G [0, 1], F(x(s)) = 0; 

(iii) for any k G {1, n T — 1}, Tk(x(s)) = for all s G [0, 1]. 

(iv) for any s G [0,1], if r k (x(s)) ^ for some k G {l,...,n T — 1}, then 
A fc+ i(«)=A^ 1 ,...,A nT ( S )=A« r 1 ; 
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Lemma 2j Let j G {l,...,n T — 1} and let F, r, U*, A*, satisfy 
the hypothesis of the Theorem. Assume that there exists an open neighbor- 
hood W® of (U*,A*) such that for any (U j , A j ) G n T, s there 
exists a continuous piecewise smooth curve x^\s) = (U^(s), A^(s)), s G 
[0,1], x ij \0) = (U j ,A j ), x®(l) = (U*,A*) with the following properties: 

(ij) x (j) ([0, 1]) C and consists of n T — j arcs, specifically, there 

exist numbers < < < ... < Sn] = 1 such that each A^ = 

x ^\[ s k-v s k D' k =j + 1, ...,7V is an arc; 
(Uj) for any s 6 [0,1], F (x^(s)) = 0; 

(Hi j) for any k G {j + 1, ...,n T }, if x G A^ , then T\(x) = ... = Tk-i(x) = 

0. 

(iVj) for any s G [0, 1], if r k (x^'(s)) ^ for some k G {1, ...,n T — 1}, 
then A k+1 {s) = A{ +1 , A nT (s) = A^. 

Then there exists an open neighborhood W^>~ 1 > of (Z7*,A*) snc/j that for 
any (W" 1 , A 3-1 ) G pt^ -1 ) n£j_i i/iere exists a continuous piecewise smooth 
curve xV-Vfa), s G [0,1], a;Cj-i) (o) = (CP -1 , A 3 ' -1 ), ^'"^(l) = (U*,A*) 
which satisfies the conditions (i) — (iv) above with j substituted by j — 1. 

Proof of Theorem. Lemma 1 guarantees the existence of a neighbor- 
hood described in the hypothesis of Lemma 2„ T _i (take W^ nT ~^ = W). The 
conclusion of the latter assures the existence of a neighborhood described 
in the assumption of Lemma 2„ T _ 2 , and so on. At the end one finds that 
there exists a neighborhood described in the conclusion of Lemma 2\. 
The properties of coincide with the requirements for V. Thus take 

V = W<®. 

Proof of Lemma 1. By the Inverse Mapping Theorem |I| [26[], applied 

to F^ 1 ^ there exist open neig hborhoods A, A'oi (U*,A*) and F^\U*,A*) = 
0, respectively, such that the restriction F^ nT ^ : A — > A' has a continuously 
differentiate inverse H : A' — > A. Let B be an open ball in Y x W lT centered 
at and contained in A'. Define W = H(B). Let (IF*- 1 , A"-" 1 ) G WnE^-i. 
Define 

x(s)= J ff( 2 /(s)) J sG[0,l], (3.14) 

where 

2 /( S ) = (y y ( S ), 2 / 1 ( S ),..., 2 /„, T ( S )) = (O y O,...,0,(l-s)r nT (^- 1 ,A"- 1 )), 

(3.15) 

and 0y is the zero element in Y. 
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To see that x satisfies condition (i), notice that y(s), s G [0,1], is a 
rectilinear segment connecting y(0) = F^ nT \U nT ^ 1 , A™ T_1 ) to y(l) = 0. The 
latter is the center of the ball B and fK)([/"^ 1 ) A nT_1 ) g £?. Therefore 
y[0, 1] C i? (so, in particular, H(y(s)), s G [0,1], is well defined). We 
conclude that x([0, 1]) = H(y([0, l]))is an arc. 

That x(s) satisfies condition fii) is clear from F(x(s)) = Fy" T ^(x(s)) = 
y Y (s)=0, se[0,l]. 

Notice that, for k — 1, ...,n r — 1, rfc(x(s)) = F^" T ^(a;(s)) = yk(s) = for 
all s G [0, 1], i.e. condition (iii) holds. 

Condition (iv) is a trivial consequence of (iii). 

Proof of Lemma 2j. By the Inverse Mapping Theorem applied to 
there exist open neighborhoods A, A', of (U*,A*) and of y* = F^(U*, A*) = 
(0y, 0, 0, A* +1 , A* r ) , respectively, such that F^) : A ^ A' has a contin- 
uously differentiable inverse H : A' -> A. Then (W^ n A) = H' 1 {W {j) n A) 
is open and contains y* (since both W {j) and A are open, both contain 
(U*,A*), and is continuous). Therefore there exists an open ball B in 
YxW 1 - centered at y* and contained in F® {W^ n A) . Set wCj-i) = #(£). 

Let us verify that W^~^ possesses the desired properties. Let A-' -1 ) G 

jyO'- 1 ) n Sj_i. Define a parametrized curve x(s) by 

z( S )=tf(y(s)), sG[0,l], (3.16) 

where 

y(s) = (y Y (s), yi (s), ...,y nT (s)) = (0y, 0, 0, (1 - s)^- 1 , A^ 1 ), A£}, A^; 1 ) 

(3.17) 

Notice that y(0) = F ( j\U j ~ 1 ,A j ~ 1 ) G B and 

II 2/(0) -y* || = max{|| yy || y , || (yi,...,y nT ) \\r^} =|| (j/i, Z/nJ ||m™t 

= O^- 1 , a- 1 )] 2 + (Aj;} - a* +1 ) 2 + ... + (ac 1 - a; t ) 2 ) 1/2 

> (( A S - a ; + i) 2 + - + (K: 1 - Kf) 1/2 =11 y(i) - v* II, (3.18) 

i.e., y(l) is not farther from the center y* of the ball B than y(0). Hence 
2/(1) G B, and so y([0, 1]), which is a rectilinear segment, is contained in B. 
This implies, first, that y([0, 1]) is in the domain of H, hence x(s) is well 
defined, and, second, that x([0, 1]) = H(y[0, 1])) C W^'- 1 ) is an arc, since H 
is a diffeomorphism and y([0, 1]) is a rectilinear segment. 
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Notice that x(l) G fl Sj. Therefore, by assumption, there exists a 
continuous piecewise smooth curve x {i) (s), s G [0, 1], such that a#)(0) = x(l), 
= (£/*, A*), and the conditions (ij) - (ivj) are satisfied. Define 

x 0'-i)f s ) - / x (2s) ? 0< S <l/2, 

X lSj ~\ xW(2s-l), 1/2 < s < 1. ^ i9j 

By construction, x^' -1 )(s) is continuous, piecewise smooth, and x^~^([0, 1]) C 
W^-V. Also ^'^([O, 1]) consists of n T - j + 1 arcs, since a^ _1 )([0, ±]) = 
x([0, 1]) is an arc, whereas a# _1) ([§, 1]) = x (j) ([0, 1]) was assumed to con- 
sist of n r — j arcs. More precisely, if s^\...,Sn} are such as assumed in 
(iiij), then we define s^~^ = 0, s^~^ = (sjj^ + l)/2, = j,...,n T (clearly, 
< < < ... < sV ] = 1). Then A^ = x^ds^ , s^]) = 

a;(j-i)([o, I]) = x([0, 1]) and, for fc = j+1, ...,n T , A^ = x^([s^\ sg _1) ]) = 

^■-^([ftkti, £l±1]) = rr^Qs^x, 4 J) 1) = ^fc } are arcs - Thus condition (i^) 
is satisfied. 

Condition (iij-i) is satisfied for x ( -^ 1 \s), since it is satisfied for both x(s) 
and x^\s). 

If x G ^4^ _1 \ i-e. x = x(s) for some s G [0, 1], then, for m — 1, j — 1, 
r m (x) = = F^(x(s)) = F^(H(y(s))) = y m {s) = 0. If x G ^'"^ 

for some A; G { j + 1, ...,n T }, then, since ^f - ^ = ^4.^, we have T\{x) = ... = 
i~k-i( x ) — by (iiij). Thus (iiij) holds. 

To verify (ivj_i) for x^ _1 ^(s), one has to consider s G [0, |) only, since, 
for s G g, 1], x^' _1 )(s) = x {j \2s - 1), and (i Vj ) applies. If s G [0, §), then 
xfr' -1 )(s) = x(2s), so, from the definition of x{s), T\ [x^~^(s)) = ... = 
Tj_i [x^~^(s)) = 0. Thus the implication in (ivj-i) has to be verified for 
k > j only. Suppose that Tj (^-^(O)) = T j (CP -1 , A J_1 ) ^ (otherwise the 
arc x([0, 1]) degenerates into a single point which has already been 

considered. Then, for any s G [0, |), 

r,(x(2 S )) = r J (H(y(2s))) = F^(H(y(2s))) = y 3 (2s) = (l-2s)r^-\ A^ 1 ) ^ 0. 

(3.20) 

Let us first verify (ivj-i) for k = j. For any x G X, its last n T — j coordinates 
coincide with those of F^\x) by the definition of c.f. (3.11a). But 
F^\x(s)) = y(s) whose last n T — j coordinates are A^\, A^" 1 , respectively. 
This proves the implication contained in (ivj_i) for k = j. 
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Validity of (ivj_i) for k > j follows immediately from the above observa- 
tion that the last n T — j coordinates of x(s), < s < 1, are A^}, A^ 1 , 
respectively. 

Remark 5 It follows from the construction of the curves x(s) in the proofs 
of Lemmas 1 and 2j that, when Aj is 'freed," the value of \ Tj(x(s)) \ strictly 
decreases, as s goes from to 1, provided Tj(x(0)) ^ 0. This is immediate 
from Tj (x{s)) = Fj j \H(y(s))) = y 3 {s) = (1 - s)r,(x(0)). 



Remark 6 Curves constructed in the proofs of Lemmas 1 and 2j have finite 
length. Indeed, they consist of a finite number of arcs of finite length, see 
Remark 4- 



4 Examples 

4.1 Example 1: Homoclinic orbits in a 3-D singular 
perturbation problem. 

We compute a homoclinic orbit for the 3-D system |Tl| 



x' = (2 - z)a(x - 2) + (z + 2) [a(x - x ) + (3(y - y )] , (4. la) 

y' = (2 -z)[d(b-a)(x-2)/4 + by] + (z + 2)[-/3(x-x ) + a(y 

z' = {A - z 2 )[z + 2 - m{x + 2)] - ecz, (4.1c) 

where a = 1, b = 1.5, c = 2, d = —.2, m = 1.1845, a = .01, f3 = 5, xo = — .1, 
y Q = —2. The parameter e is taken as variable. In this case n = 2 and ni = 1 
in Equation (2.1). The discretization is orthogonal collocation with piecewise 
polynomials, using 25 subintervals and 4 collocation points per interval. A 
relative Newton tolerance of 10 -8 is used for u and A. 

Step 0. Initialize the problem parameter: e = .01, and the algorithm 
parameters: e = 10 -4 , T = 10~ 2 , C\ = 1, c 2 = 0. 

Step 1. Compute a branch of solutions to the equations (2.4), (2.7), 
(2.11), (2.12a), and 

d o = J2c jqoj , (4.2) 
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with scalar variables T, E\ G IR, do, d\ G IR 3 , rx,r 2 G R, in the direction 
of increasing T until T\ crosses zero. There are ten scalar variables in this 
continuation. Final values are T = 2.7568, e\ = 3.509, and r 2 = .0223. 

Step 2. Fix n = and free c\ and c 2 in (4.3). Compute a solution 
branch to the system (2.4), (2.7), (2.11), (2.12), and (4.3), in the direction of 
decreasing c 2 until r 2 crosses zero. There are now eleven free scalar variables. 
Final values are c x = 1.0, c 2 = -1.51 x 10~ 4 , T = 2.7580, and e t = 2.548. 

Step 3. The free parameters ci, c 2 G K are replaced by do G IR 3 ; see Equa- 
tion (2.6). Fix T\ — and free e G K and hence the matrices Qo, Qi, To, T\ G 
IR 3 x IR 3 . Compute a solution branch to (2.4), (2.5a), (2.6), (2.7), (2.9) - 
(2.11) in the direction of decreasing e% (to increase the accuracy of the or- 
bit) until €i = to- There are 44 free scalar variables. Terminal values are 
e = .009333142, T = 2.7679, and e x = 10" 4 . 

Note, that it took only 2 steps to get n = r 2 = 0, while the general 
algorithm requires 3 steps to accomplish this. In particular, since T\ crossed 
zero already in the first step, we were able (without changing the total number 
of free scalar variables) to keep T free in the second step, while T\ — was 
fixed. 



4.2 Example 2: Heteroclinic orbits in a 3-D Josephson 
Junction problem. 

A singularly perturbed sine-Gordon equation, modeling magnetic flux quanta 
( "fluxons" ) in long Josephson tunnel junctions with nonzero surface impedance, 

M'\i) - (1 - c 2 )0"(O - ac0'(O + sin 0(0 -7 = 0, (4.3) 



was studied by several authors, see e.g. || and references therein. In [14 
we computed single and multiple fluxon solutions, which are heteroclinic 
orbits (or homoclinic orbits on a cylinder). The algorithm in [14] does not 



allow computation of the transition from two real eigenvalues to a complex 
conjugate pair. Here we accurately compute this transition point. The three- 
dimensional first order system is 

0'i = 2 , (4.4a) 
02 = 03, (4-4b) 
<j)' 3 = [(l-c 2 )0 3 + CH^-sin^ + i\//3c (4.4c) 
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We compute a heteroclinic orbit with u = (arcsiri7, 0, 0)and U\ = (arcsin7 + 
7r, 0, 0). Throughout, a = .18 and (3 = .1 are kept fixed, and 7 and c vary. 
In this case n = 1 and n\ = 2 in Equation (2.1). Discretization is as in the 
preceding example. 

Step 0. Initialize the problem parameters, 7 = .1, c = .6, and the 
algorithm parameters, e = 10~ 4 , e\ = .6283, T = 10" 2 , and c\ = 1. 

Step 1. Compute a solution branch to (2.4), (2.7), (2.11), (2.12a), 
(4.2), with scalar variables T, e\ G R, do,d\ G M 3 ,r G M, in the direction 
of increasing T until r crosses zero. There are nine free scalar variables. 
Terminal values are T = 9.336, and e 1 = 3.642. 

Step 2. Fix Ti = 0. Free c G M, and hence the matrices Q , Qi G M 3 x M 3 , 

and the entries 1? 2 ) ^?3' ^22; ^23; ^33 an d ^i2> ^13' ^22' ^23' ^33 °f the matrices 
Tq,T\ G 1R 3 x M 3 , respectively, as well as tifa, in Equation (4.5) below. 
Compute a branch of solutions to the system (2.4), (2.5a), (2.6) - (2.10), and 



in the direction of decreasing ei, until e\ — e . This step involves 44 free scalar 
variables. Terminal values are c = .3404, T = 21.13, t\ = 10~ 4 , /% = 1.073, 
/i n = —2.600, and /i 12 = —1.047. Above, t® 2 = Re/i 02 , t® 3 = Re/i 03 , = 
Re/x 11 , t 2 2 = R- e / i i2- The parameters ti,t 2 in (4.5) serve as test functions that 
cross zero when /i 02 and /i 03 (and /i u and /i 12 ) are multiple. 

Step 3. Fix ei and free 7 G M. Compute a solution branch with pa- 
rameters 7 and c in the direction of decreasing 7 until ti and t 2 cross zero. 
The equations are as in Step 2. Terminal values are 7 = 4790, c = .3404, 
T = 12.92, /i 01 = 1.622, and /jl u = /jl 12 = —2.534. The double eigenvalue was 
located with accuracy 10~ 8 . 

Step 4. To continue the complex conjugate pair of the eigenvalues, fix 
ti = t 2 = (Re/x n = Re/x 12 ) and free ^32, ^21 (I m / i o2J m / i i2)- Compute a 
solution branch with parameters 7 and c to the same system as in Step 3 in 
the direction of decreasing 7. Final values are 7 = 0.8830, c = 1.000, and 
T = 30.67. 

Note, that it took only one step to get r = 0, while the general algorithm 
requires two steps to accomplish this. 






33' 
1 



(4.5a) 
(4.5b) 
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4.3 Example 3: Heteroclinic orbits in a modified 3-D 
Josephson Junction problem. 

Consider system (4.4) with reversed time, 

01 = -02, (4.6a) 

02 = -03, (4.6b) 

3 = -[(1-c 2 ) 03 + ac0 2 -sin0 1 +7] //5c. (4.6c) 

We compute a heteroclinic orbit with u = (arcsin7 + ir, 0, 0) and -Ux = 
(arcshi7, 0, 0). As before, a = .18 and (3 = .1 are kept fixed, and 7 and 
c vary. The problem parameters are chosen so that n = 2 and ni = 1 in 
Equation (2.1), where /% and /% are a complex conjugate pair of eigenvalues. 
Discretization is as in the preceding example. 

In this problem the starting direction of trajectories near Uq is unknown, 
making it difficult to generate a starting orbit in Step with u(l) in a small 
neighborhood of u±. Indeed, this example also illustrates the more global 
applicability of the successive continuation approach. 

Step 0. Initialize the problem parameters, 7 = .608, c = —.95, and the 
algorithm parameters, e = 10~ 4 , e\ = .6283, T = 10~ 2 , and c\ = 1, c 2 = 0. 
Initially Re/% = Re/% = 1.508, Im/% = — Im/% = 1.388. 

Step 1. Compute a solution branch to (2.4), (2.7), (2.11), (2.12a), 
(4.2), with scalar variables T, e\ GR, do, d\ G M 3 , 77, r 2 G K, in the direction 
of increasing T until 7~2 crosses zero. There are ten free scalar variables. 
Terminal values are T = 6.098, and e x = 6.1738, 77 = —.9500. 

Step 2. Compute a solution branch to (2.4), (2.7), (2.11), (2.12), (4.2), 
with scalar variables T,e 1 G K, c% ^1 G ffi 3 , Ci,c 2 ,Ti G M, in the direction 
of increasing T. There are nine free scalar variables. Terminal values are 
T = 9.069, ei = 1.737, and 77 = -.3718. 

Step 3. Fix ei. Free cGR, and hence the matrices Qo,Qi G M 3 x M 3 , 
and the entries t^i, ^12, ^13, ^21, 4, 4, 4 an d t u , t 12 , t l3 , t 2 2, ^23, ^32, ^33 °f the 
matrices T ,Ti G M 3 x IR 3 , respectively. Compute a branch of solutions to 
the system (2.4), (2.5a), (2.6) - (2.10), and 

tl = ^-4 = 0, (4.7a) 

t 2 = 4-4 = 0, (4.7b) 

in the direction of decreasing |ri| until 77 crosses zero. This step involves 
44 free scalar variables. Terminal values are c = —.8995, T = 7.436. Above, 

4 = R e ^01, 4 = R e A t 02, 4 = R e / i 12, 4 = R e A t 13- 
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Step 4. Fix T\ = and free 6\ G HL Compute a branch of solutions in 
the direction of decreasing ex, until e x = e .The equations are as in Step 3. 
Terminal values are c = —.9027, T = 13.16, e\ = 10~ 4 . 

Note that, as in the general algorithm, it took three steps to get T\ = 

T 2 = 0. 

4.4 Example 4: Heteroclinic orbits in a 4-D singular 
perturbation problem. 

The existence of traveling wave front solutions to the singularly perturbed 
react ion- diffusion system 

v t = v xx + v(v - a)(l - v) - w, w t = 5w xx + e(v - 710), (4.8) 



for small positive e and 5, was established by Deng ||12|| . In moving coordi- 
nates, Vx = v(z), v 2 = v (z), Wx = w{z), w 2 = w'(z) with z = t + cx, the 
reduced ODE is 

v[ = v 2 , (4.9a) 

v' 2 = cv 2 — vx(l — vx)(vx — a) +wx, (4.9b) 

w[ = w 2 , (4.9c) 

w 2 = [cw 2 — e(vx — jWx)]/6. (4.9d) 

We compute a heteroclinic orbit with u = (0, 0, 0, 0) and ux — (| + 3° + 
i v /(7-2a7 + 7a 2 -4)/7, 0, Vxh, 0), with 5 = e = .001 and a = .3 fixed 
and and 7 and c variable. Initially 7 = 13.8 and c = .257. In this case 
ux = (.8736878,0,-0633107,0), and n = n x = 2, where the relevant eigen- 
values are fi 01 =.6957, /i 02 = 257.0537, \i x x = —.4415, and /i x 2 = —0.0668. 
Throughout we use a discretization with 50 subintervals, 4 collocation points 
per interval and relative Newton tolerances 10 -8 for u and A. 

In this problem there is strong divergence of trajectories near uq, making 
it difficult to generate a starting orbit in Step with u(l) in a small neighbor- 
hood of Ux- Indeed, this example also illustrates the continued applicability 
of the algorithm to such cases. 

Step 0. Initialize the problem parameters, 7 = 13.8, c = .2570, and the 
algorithm parameters, e = .6, ex = .5161, T = 10" 5 , Ci = 1, c 2 = 0. 

Step 1. Compute a solution branch to (2.4), (2.7), (2.11), (2.12a), 
(4.3), with scalar variables T,e 1 ,r 1 ,r 2 G M, d ,dx G 1R 4 , in the direction 
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of increasing T until T\ crosses zero. Terminal values are T = .0472 and 
ei = .5056. 

Step 2. Fix T\ = 0, free c\ and c 2 , and compute a solution branch 
to (2.4), (2.7), (2.11), (2.12), (4.3), in the direction of increasing T, with 
scalar variables T, ei, r 2 , ci, c 2 G M, di G M 4 . Final values are T = 1.363, 
e x = .4024, ci = 1.00000, c 2 = 1.9 x 10" 8 , and r 2 = .521. 

Step 3. The scalar variables and equations are as in Step 2, except that 
T is fixed and eo is free. Compute a branch of solutions in the direction of 
decreasing eo and t\ to locate zero of r 2 . Terminal values are eo = .3199, t\ = 
.3991, c\ = 1.00000, and c 2 = 1.3 x 10~ 8 . There is a very sensitive dependence 
on the "shooting angle", represented by c 2 . 

Step 4. The free parameters Ci,c 2 G K. are replaced by d G M 4 . Fix 
r 2 = and free T. Also free the problem parameter c G M, and hence 
the matrices Qo,Qi,Tq,Ti G M 4 x R 4 . Compute a solution branch to the 
system (2.4)-(2.10) in the direction of increasing T. Terminal values are 
c = .3437209, T = 12.66, and e = 10~ 4 . 

Step 5. The scalar variables and equations are as in Step 4, except that 
eo is fixed and e x is free. Compute a solution branch in the direction of in- 
creasing T. Final values are c = .2572501, T = 71.31, ei = 5 x 10~ 3 . 



5 Discussion 

The initial detection of a connecting orbit in a dynamical system is generally 
a difficult task; often with extremely sensitive dependence on initial condi- 
tions and on parameter values. This is even more so in the case of singularly 
perturbed equations. Thus there is use for a systematic procedure that has 
a good chance of success, even in difficult problems. In this paper we have 
presented a successive continuation method for locating connecting orbits. 
We have shown that the procedure works, provided the connecting orbit is 
isolated and the initial orbit sufficiently close. Thus our analysis only guar- 
antees local convergence of the method, even though it has been designed 
for extended convergence properties, as is well illustrated by our numerical 
examples, which include some "hard" problems. A complete presentation of 
global convergence properties is beyond the scope of this paper, and perhaps 
best presented in a more general context. The key ideas in a global analysis 
of the successive continuation algorithm are the following. It is assumed that 
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the problem can be reduced to that of finding a small number of parameters 
for which an equal number of equations is to be satisfied. Given a sub- 
set of the equations that are satisfied, a one-dimensional continuum of such 
solutions is followed, and points where any of the remaining equations are 
satisfied are accurately located. Proceeding, inductively, one then continues 
an enlarged set of satisfied equations. At any stage of the algorithm, the 
continuation procedure works locally, provided the solution points are "reg- 
ular" ]TJJ. Generically this is the case, while for nongeneric problems (which 
are often encountered!) one can regularize the problem by adding (if neces- 
sary) unfolding parameters and (again, if necessary) regularizing equations. 
In fact, the method of pseudo-arclength continuation itself is the simplest 
nontrivial example of this procedure. 

The power of the successive continuation procedure is then its ability to 
reach a solution from a far away starting point, and, in fact, to locate multiple 
solutions and often all solutions, provided the regularity assumptions are 
satisfied and provided the solution(s) are reachable. It is the latter condition 
that need not always be satisfied. In fact, it is easy to construct simple 
(algebraic) examples where a solution is not reachable from a given starting 
point. One could argue that this problem can also be solved by adding 
unfolding parameters, but in practice it is not often clear how to do this, for 
example, in the case of computing connecting orbits. Nevertheless, as our 
numerical experience has shown, including the examples presented in this 
paper, the successive continuation method provides, at least, a useful tool 
that often gives results where other methods fail. 

Note that, in our computation of connecting orbits, the above reduction 
to a "small-dimensional problem" does not require the problem to be posed as 
a "shooting problem" , which would make the algorithm useless, for example, 
for the singular perturbed equation in Section 4.1. Throughout, the small 
dimensional problem remains embedded in the full set of equations, which 
are solved by continuation in the full space. Even the initial "integration" 
(see Algorithm, Step 1) is done by continuation of complete orbits, in order 
to be able to monitor, and react to (e.g. adaptive mesh refinement), sensitive 
dependence, for example, on the left boundary conditions. 

Although, as pointed out above, a discussion of global convergence prop- 
erties is perhaps best presented in terms of the global topology of the under- 
lying manifolds, one can also obtain global results by making certain global 
assumptions on the vector field, although this is practically less useful in 
the present context. In fact, applying the theory in pB| , one can show the 
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following. Let Ff^ T ' 
of :X^Yx\ 



X h -+Y h x R n - be the discretization |§ pp. 748-749] 
1 Ut defined by (3.11b). Let Q C X^ be a bounded open 
set with a smooth connected boundary dfl. Suppose that the operator F^ 1 



satisfies Smale boundary conditions |25[ 



a) DuFfr\Uh, Ah) is nonsingular on dQ; and either 



b) {DuFtHUH^))- 1 ^' 
^(DuF^iU^An))-^' 
Then by a slight generalization of [25, Th. 2.4] one can show that V£/£ G 
dfl there exists a piecewise smooth path with starting point (?7°,A°) and 



^Ft\U h ,A h) 
^F^\u h ,A h ) 



points into Q, V£4 G <9f2; or 
points out of Q, \/Uh G dfl. 



terminal point (U^,A* h ) with F^' (U£, A* h ) = 0. 
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